setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210321.H3K分析/")
library(ggplot2)
library(ggsci)

group1=c("117012","13649","8544","807","200")
group3=c("enhancer","promoter","all","")
i=1
j=1
k=1
fn=paste0(group1[i],"/","consider.all.VAF.Brain_Hippocampus_Middle.",group3[k],"_statis.csv")
tmp=read.csv(fn,header = T)
group2=unique(as.character(tmp$X))

p <- list()
for (j in 1:length(group2)) {
  rt <-tmp[-1,]
  for (i in 1:length(group1)) {
    for (k in 1:length(group3)) {
      fn=paste0(group1[i],"/","consider.all.VAF.Brain_Hippocampus_Middle.",group3[k],"_statis.csv")
      tmp=read.csv(fn,header = T)
      tmp=tmp[tmp$X==group2[j],]
      tmp$Term=group1[i]
      tmp$H3k=group2[j]
      tmp$regions=group3[k]
      rt=rbind(rt,tmp)
      }
  }
  rt.plot=rt
  rt.plot[rt.plot$regions=="",]$regions = "Not enhancer or promoter"
  rt.plot$oppo.ratio=rt.plot$opposite/(rt.plot$same+rt.plot$opposite)
  rt.plot$Group <- paste0(rt.plot$Term, "_",rt.plot$H3k)
  rt.plot$Sign <- ifelse(rt.plot$pvalue<0.05,"S","N.S")
  rt.plot=data.frame(Same=rt.plot$same_ratio,Opposite=rt.plot$oppo.ratio,Group=rt.plot$Group,Sign=rt.plot$Sign,regions=rt.plot$regions,H3k=rt.plot$H3k)
  rt.plot1=rt.plot[,-1]
  colnames(rt.plot1)[1] <- "ratio"
  rt.plot1$Ratio <- "Opposite"
  
  rt.plot2 <- rt.plot[,-2]
  colnames(rt.plot2)[1] <- "ratio"
  rt.plot2$Ratio <- "Same"
  
  rt.plot <- rbind(rt.plot1, rt.plot2)
  rt.plot$regions <- factor(rt.plot$regions, levels = c("enhancer", "promoter", "Not enhancer or promoter", "all"))
  rt.plot$Sign <- factor(rt.plot$Sign, levels = c("S", "N.S"))
  rt.plot$Ratio <- factor(rt.plot$Ratio, levels = c("Opposite", "Same"))
  rt.plot$Group <- factor(rt.plot$Group, levels = paste0(group1,"_",group2[j]))
  p[[j]]=ggplot(data = rt.plot, 
         aes(x = Ratio, y = ratio, group = regions, color = regions))+
    geom_line(size=1)+geom_point(aes(shape=Sign),size=3)+
    #scale_shape_manual(values = c(1,2)) +
    theme_light()+labs(x="",y="Ratio",title = "")+
    theme(axis.title.y=element_text(size=16),axis.text.y=element_text(face="bold",size=15),axis.text.x=element_text(size=15))+
    scale_color_npg() +
    facet_wrap(~Group, scales = "free")

}
